Final study sample
participants <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/sample_info/7T_MP2RAGE_finalsample_demographics.csv") #215 scans from 140 subjects. csv generated by /sample_construction/finalsample_7Tmyelin.RmdGlasser regions list
glasser.regions <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist.csv")
glasser.snr.exclude <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/SNR/glasser_SNR_exclusion.csv")
glasser.frontal <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist_frontallobe.csv")glasser.plotting <- glasser.regions
glasser.plotting$cortex <- "cortex"
glasser.plotting$cortex[glasser.plotting$orig_parcelname %in% glasser.snr.exclude$orig_parcelname] <- "exclude"
glasser.plotting <- glasser.plotting %>% select(label, cortex)
glasser.plotting <- rbind(glasser.plotting, data.frame(label = "rh_???", cortex = "exclude"))
glasser.plotting <- rbind(glasser.plotting, data.frame(label = "lh_???", cortex = "exclude"))
glasser.plotting <- glasser.plotting %>% filter(label != "lh_L_10pp")Depth-dependent R1 in HCP-MMP regions
#R1 by region matrices at 7 cortical depths
myelin.glasser.7T <- readRDS("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/depthR1_glasseratlas_finalsample.RDS") #generated by /surface_metrics/surface_measures/extract_depthdependent_R1.Rordered_depths <- c("depth_7", "depth_6", "depth_5", "depth_4", "depth_3", "depth_2", "depth_1")
depth_colorbar <- c("#004A38FF", "#006B63FF", "#008FA7FF", "#6992CC", "#A29DC4", "#C0BFE3", "#E2CFE5")Regional average R1 in HCP-MMP regions
#Calculate across-depths average regional R1
regionalaverage.myelin <- do.call(rbind, myelin.glasser.7T) #R1 at all depths
regionalaverage.myelin <- regionalaverage.myelin %>% select(contains("ROI"))
regionalaverage.myelin <- colMeans(regionalaverage.myelin) %>% as.data.frame() %>% set_names("mean.R1")
regionalaverage.myelin <- regionalaverage.myelin %>% mutate(orig_parcelname = row.names(regionalaverage.myelin))Myelin Maps
#Myelin Water Fraction Imaging
MWF.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/MyelinWaterFraction_Liu2019/source-Liu2019_desc-MWFmap_space-fsaverage_den-164k_atlas-glasser360.pscalar.nii")
MWF <- data.frame(MWF = MWF.cifti$data, orig_parcelname = names(MWF.cifti$Parcel))
#Myelin Basic Protein Gene Expression
MBP.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/AHBA_magicc/expression_maps/source-magicc_desc-MBPexpression_space-fsLR_den-32k.pscalar.nii")
MBP <- data.frame(MBP = MBP.cifti$data, orig_parcelname = names(MBP.cifti$Parcel))
#Sensorimotor-Association Axis
SAaxis.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii")
SAaxis <- data.frame(SA.axis = rank(SAaxis.cifti$data), orig_parcelname = names(SAaxis.cifti$Parcel))df.list <- list(regionalaverage.myelin, MWF, MBP, SAaxis) #dfs to merge
myelin.maps <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list)
myelin.maps <- merge(myelin.maps, glasser.regions, by = c("orig_parcelname"), sort = F)
myelin.maps <- myelin.maps[!(myelin.maps$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels
write.csv(myelin.maps, "/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_corticalmyelin_anatomy/myelin_maps/myelin_maps.csv", quote = F, row.names = F)myelin.colorbar <- c("#edf6ff", "#c2d4ed", "#809dc4", "#345c91", "#0e4669")myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = mean.R1), , colour=I("gray50"), size=I(.06)) +
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.54, 0.58), oob = squish, na.value = "gray75") + theme(legend.position = "none")ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/R1_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)MWF Map
myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = MWF), colour=I("gray50"), size=I(.06)) +
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.02, 0.05), oob = squish, na.value = "gray75") + theme(legend.position = "none")ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MWF_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)Correlation
cor.test(myelin.maps$mean.R1, myelin.maps$MWF, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: myelin.maps$mean.R1 and myelin.maps$MWF
## S = 2261576, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6422757
Correlation Plot
myelin.maps %>% filter(MWF < 0.08) %>% #remove 2 MWF major outliers in PHA1
ggplot(., aes(x = MWF, y = mean.R1)) +
geom_point(color = c("#8AABD5FF"), size = 0.2) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
labs(x="\nMyelin Water Fraction", y="R1\n") +
scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
theme(
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MWF_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)MBP Map
myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = MBP), colour=I("gray50"), size=I(.06)) +
scale_fill_gradientn(colors = myelin.colorbar, limits = c(-0.65, 1), oob = squish, na.value = "gray75") + theme(legend.position = "none")ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MBP_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)Correlation
cor.test(myelin.maps$mean.R1, myelin.maps$MBP, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: myelin.maps$mean.R1 and myelin.maps$MBP
## S = 379704, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5279897
Correlation Plot
myelin.maps %>%
ggplot(., aes(x = MBP, y = mean.R1)) +
geom_point(color = c("#8AABD5FF"), size = 0.2) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
labs(x="\nMyelin Basic Protein Expression (z)", y="R1\n") +
scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
theme(
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MBP_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)S-A axis
myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = SA.axis), colour=I("gray50"), size=I(.06)) +
scale_fill_gradientn(colors = myelin.colorbar, trans = 'reverse', na.value = "gray75") + theme(legend.position = "none")ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/SAaxis_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)Correlation
cor.test(myelin.maps$mean.R1, myelin.maps$SA.axis, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: myelin.maps$mean.R1 and myelin.maps$SA.axis
## S = 10280562, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.6261257
Correlation Plot
myelin.maps %>%
ggplot(., aes(x = SA.axis, y = mean.R1)) +
geom_point(color = c("#8AABD5FF"), size = 0.2) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
labs(x="\nSensorimotor-Association Axis", y="R1\n") +
scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
theme(
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/SAaxis_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)#Compute average frontal cortex volume fraction at each depth for final study participants
cortexpve <- read_parquet("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/7T_BrainMechanisms_cortexpve.parquet") #parquet with cortex fraction measures
cortexpve <- cortexpve %>% filter(atlas == "glasser") #just for glasser atlas
cortexpve <- cortexpve[!(cortexpve$StructName %in% glasser.snr.exclude$orig_parcelname),] #just for good regions
cortexpve <- cortexpve[cortexpve$StructName %in% glasser.frontal$orig_parcelname,] #just for frontal lobe
extract_depth_cortexpve <- function(mymeasure){
measure.df <- cortexpve %>% select(subject_id, session_id, StructName, all_of(mymeasure))
names(measure.df) <- c("subject_id", "session_id", "orig_parcelname", "cortex.fraction")
depth.cortexpve <- merge(measure.df, participants, by = c("subject_id", "session_id")) %>% select(cortex.fraction) %>% colMeans(.$cortex.fraction) #get cortex.fraction for final study participants only, and then average for this depth
return(depth.cortexpve)}
pve.measures <- list("Mean_cortexpve.1.00%", "Mean_cortexpve.0.9%", "Mean_cortexpve.0.8%","Mean_cortexpve.0.7%","Mean_cortexpve.0.6%","Mean_cortexpve.0.5%","Mean_cortexpve.0.4%","Mean_cortexpve.0.3%","Mean_cortexpve.0.2%","Mean_cortexpve.0.1%", "Mean_cortexpve.0.0%") #measures to extract data for
cortexpve.frontallobe <- lapply(pve.measures, function(x) {
extract_depth_cortexpve(x)})
cortexpve.frontallobe <- do.call(rbind, cortexpve.frontallobe) %>% as.data.frame()
cortexpve.frontallobe$depth <- c("a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k")
cortexpve.frontallobe$depth <- factor(cortexpve.frontallobe$depth, ordered = T, levels = c("k", "j", "i", "h", "g", "f", "e", "d", "c", "b", "a"))
kable(cortexpve.frontallobe)| cortex.fraction | depth |
|---|---|
| 0.7103157 | a |
| 0.8388772 | b |
| 0.9281964 | c |
| 0.9749537 | d |
| 0.9922900 | e |
| 0.9960812 | f |
| 0.9925421 | g |
| 0.9710414 | h |
| 0.9071242 | i |
| 0.7684002 | j |
| 0.5157232 | k |
cortexpve.frontallobe %>%
ggplot(aes(x = cortex.fraction, y = depth)) +
geom_point(color = "#809dc4", shape = 15, size = 2) +
theme_classic() +
xlab("\nCortex Fraction") +
ylab("Cortical Depth\n") +
theme(
legend.position = "none",
axis.text.x = element_text(size = 6, family = "Arial", color = c("black")),
axis.text.y = element_blank(),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks.x = element_line(linewidth = .2),
axis.ticks.y = element_blank()) ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1_supplementary/Depthplot_cortexfraction_frontal.pdf", device = "pdf", dpi = 500, width = 2.5, height = 2.5)#Compute mean R1 in each region at each cortical depth
regional_means <- function(input.df){
meanmap <- input.df %>% select(contains("ROI")) %>% colMeans() %>% as.data.frame() %>% set_names("mean.R1")
meanmap$orig_parcelname <- row.names(meanmap)
meanmap <- merge(meanmap, glasser.regions, by = c("orig_parcelname"), sort = F)
return(meanmap)
}
regional.myelin.alldepths <- lapply(myelin.glasser.7T, function(depth) {
regional_means(depth)
})regional.myelin.alldepths <- do.call(rbind, regional.myelin.alldepths) #regional mean R1 at all depths
regional.myelin.alldepths$depth <- substr(row.names(regional.myelin.alldepths), 1, 7) #assign depth variable
regional.myelin.alldepths$depth <- factor(regional.myelin.alldepths$depth, ordered = T, levels = ordered_depths)
regional.myelin.alldepths <- regional.myelin.alldepths[!(regional.myelin.alldepths$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels
regional.myelin.alldepths <- regional.myelin.alldepths[regional.myelin.alldepths$orig_parcelname %in% glasser.frontal$orig_parcelname,] #only studying frontal
regional.myelin.alldepths %>%
ggplot(aes(x = mean.R1, y = depth)) +
geom_violin(aes(fill = depth), color = "white") +
stat_summary(fun = median, geom = "point", shape = 23, size = 2.25, color = "white") +
theme_classic() +
scale_fill_manual(values = depth_colorbar) +
scale_x_continuous(limits = c(0.49, 0.65), expand = c(0,0)) +
scale_y_discrete(expand = c(0, 0)) +
xlab("\nR1") +
ylab("Cortical Depth\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Depthplot_regionalR1_frontal.pdf", device = "pdf", dpi = 500, width = 2.6, height = 4.9)for(this.depth in c(unique(regional.myelin.alldepths$depth))){
depth.data <- regional.myelin.alldepths %>% filter(depth == this.depth) %>% filter(label != "lh_L_10pp")
depth.plot <- ggplot() +
geom_brain(data = depth.data, atlas = glasser, mapping = aes(fill = mean.R1), colour=I("gray10"), size=I(.08)) + theme_classic() + theme(legend.position = "none") +
paletteer::scale_fill_paletteer_c("grDevices::PuBuGn", direction = -1, limits = c(0.48, 0.61), oob = squish, na.value = "white") + theme(legend.position = "none") +
theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
new_scale_fill() +
geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex), colour=I(alpha("gray10", 0)), size=I(0)) +
scale_fill_manual(values = c(alpha("white", 0), "gray75"))
print(depth.plot)
ggsave(filename = sprintf("/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/%s_corticalmap_frontal.pdf", this.depth), device = "pdf", dpi = 500, width = 4.95, height = 2.25)
}Group Average Profiles
examplar.areas.alldepth <- rbind(area4.myelin.alldepths, area46.myelin.alldepths, areaa24.myelin.alldepths)
ggplot() +
geom_smooth(data = examplar.areas.alldepth, aes(x = mean.R1, y = depth, group = orig_parcelname, color= orig_parcelname), se = F, linewidth = .8) +
scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
theme_classic() +
xlab("\nR1") +
ylab("Cortical Depth\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Area_profiles_group.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)Individual Profiles
#R1 values in all regions at all depths for all participants
myelin.glasser.7T.alldepths <- do.call(rbind, myelin.glasser.7T)
myelin.glasser.7T.alldepths$depth <- substr(row.names(myelin.glasser.7T.alldepths), 1, 7) #assign depth variable
myelin.glasser.7T.alldepths$depth <- factor(myelin.glasser.7T.alldepths$depth, ordered = T, levels = ordered_depths)
#Wide to long format
cols_to_pivot <- names(myelin.glasser.7T.alldepths)[grepl("ROI", names(myelin.glasser.7T.alldepths))]
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% pivot_longer(cols = all_of(cols_to_pivot), names_to = "orig_parcelname", values_to = "R1")
#Select exemplar regions of interest
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI" | orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI" | orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI")
#Format grouping variables
myelin.glasser.7T.alldepths$subject_id <- as.factor(myelin.glasser.7T.alldepths$subject_id)
myelin.glasser.7T.alldepths$orig_parcelname <- as.factor(myelin.glasser.7T.alldepths$orig_parcelname)
#Plot for the oldest 3 participants
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(age > 32)ggplot() +
geom_smooth(data = myelin.glasser.7T.alldepths, aes(x = R1, y = depth, color = orig_parcelname, group = interaction(orig_parcelname, subject_id)), se = F, linewidth = .4) +
scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
theme_classic() +
xlab("\nR1") +
ylab("Cortical Depth\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Area_profiles_individual.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)regional.myelin.alldepths.wide <- regional.myelin.alldepths %>% pivot_wider(id_cols = orig_parcelname, names_from = depth, values_from = mean.R1) %>% select(-orig_parcelname)
regional.myelin.alldepths.wide <- regional.myelin.alldepths.wide %>% select(depth_7, depth_6, depth_5, depth_4, depth_3, depth_2, depth_1)
depth.R1.cors <- cor(regional.myelin.alldepths.wide) # compute correlation matrix
ggcorrplot(corr = depth.R1.cors, method = "square", type = "upper", show.diag = T, lab = F, lab_size = ) +
scale_fill_gradientn(colors = c("#edf6ff", "#c2d4ed", "#97adcc", "#345c91"))ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Depth_correlationmatrix_frontal.pdf", device = "pdf", dpi = 500, width = 4.4, height = 2.2)